% 2015 PLAGIOCLASE-LIQUID HYGROMETER
%WATERS, L.E. AND LANGE, R.A. for American Mineralogist
 
%READ ME
% The following code contains the equations to calculate the activities of
% anorthite and albite in plagioclase, dissolved wt% H2O in melt, and can be
% used to calculate pre-eruptive temperature.
 
% The spreadsheet cannot be configured to solve for temperature directly,
% because the activities of anorthite and albite are configured as a
% function of temperature, therefore, calculations of temperatures must be
% iterative.
 
% TO USE: Please crease a .csv or .xlsx file entitled "input_data" and save
% it to the directory in which this code is saved. The input_data file must
% contain the following 13 values, in the following order, with no headers:

% (1)  XAn
% (2)  XAb
% (3)  TempC (in Celsius)
% (4)  Pbars
% (5)  wtSiO2
% (6)  wtTiO2
% (7) wtAl2O3
% (8) wtFe2O3T 
% (9) wtFeOT
% (10) wtMgO
% (11) wtCaO
% (12) wtNa2O
% (13) wtK2O


% Please also create a second file entitled "results.csv". This is where 
% the results will be reported, without headers, in the following order:

% (1)  XAn,     (2)  XAb,     (3)  Aan,    (4) Aab,     (5) Temp C,
% (6)  Pbars,   (7)  wtSiO2,  (8)  wtTiO2, (9) wtAl2O3, (10) wtFe2O3T, 
% (11) wtFeOT,  (12) wtMgO,   (13) wtCaO,  (14) wtNa2O, (15) wtK2O, 
% (16) wtH2Ocalc
 
% Please note that NEW VALUES returned by the MatLab Code are:
 
% Aan (the activity of the anorthite component in plagioclase) in Column 3,
% Aab (the activity of the albite component in plagioclase) in Column 4,
% wtH2Ocalc (the calculated wt% H2O in equilibrium with the melt in Column
% 16).
 
% Please do not hesitate to contact Laura Waters if any assistance is
% needed at the following email address (laura.e.waters@gmail.com).
% Rebecca Lange may also be reached at becky@umich.edu.
 
% Thanks for your interest!





% 2015 PLAGIOCLASE-LIQUID HYGROMETER MATLAB CODE

%Clear all saved variables
clear all;

% Load data
input_data = csvread('input_data.csv');
%input_data = xlsread('input_data.xlsx');

% Define variables
XAn      = input_data(:,1);  %  molfraction An in plagioclase, 50 mol% An= 0.50
XAb      = input_data(:,2);  %  molfraction Ab in plagioclase, 35 mol% An= 0.35
TempC    = input_data(:,3);  %  Temperature Input
Pbars    = input_data(:,4);  %  Pressure input in bars, if unknown try 1000 bars
wtSiO2   = input_data(:,5);  %  Input wt% SiO2 of glass or whole rock composition
wtTiO2   = input_data(:,6);  %  Input wt% TiO2 of glass or whole rock composition
wtAl2O3  = input_data(:,7);  %  Input wt% Al2O3 of glass or whole rock composition
wtFe2O3T = input_data(:,8);  %  Input wt% Fe2O3T of glass or whole rock composition
wtFeOT   = input_data(:,9);  %  Input wt% FeOT of glass or whole rock composition
wtMgO    = input_data(:,10); %  Input wt% MgO of glass or whole rock composition
wtCaO    = input_data(:,11); %  Input wt% CaO of glass or whole rock composition
wtNa2O   = input_data(:,12); %  Input wt% Na2O of glass or whole rock composition
wtK2O    = input_data(:,13); %  Input wt% K2O of glass or whole rock composition
wtH2O    = input_data(:,14); %  Input wt% H2O of glass or whole rock composition if known

%If you have no data for any oxide component please put a "0" in its place.


% Calculate moles
molSiO2  =  wtSiO2   ./  60.085; 
molTiO2  =  wtTiO2   ./  79.899;
molAl2O3 =  wtAl2O3  ./ 101.961;
molFe2O3 =  wtFe2O3T ./ 159.692;
molFeO   =  wtFeOT   ./  71.846;
molMgO   =  wtMgO    ./  40.304;
molCaO   =  wtCaO    ./  56.079;
molNa2O  =  wtNa2O   ./  61.979;
molK2O   =  wtK2O    ./  94.195;
molAl2O3v2 = molAl2O3 - molK2O;
summoles = molSiO2 + molTiO2 + molAl2O3v2 + molFe2O3 + molFeO + molMgO + molCaO + molNa2O + molK2O;

%fprintf('%s %12.6f\n','molSiO2',molSiO2(1))
%fprintf('%s %12.6f\n','molTiO2',molTiO2(1))
%fprintf('%s %12.6f\n','molAl2O3',molAl2O3(1))
%fprintf('%s %12.6f\n','molFe2O3',molFe2O3(1))
%fprintf('%s %12.6f\n','molFeO',molFeO(1))
%fprintf('%s %12.6f\n','molMgO',molMgO(1))
%fprintf('%s %12.6f\n','molCaO',molCaO(1))
%fprintf('%s %12.6f\n','molNa2O',molNa2O(1))
%fprintf('%s %12.6f\n','molK2O',molK2O(1))
%fprintf('%s %12.6f\n','molAl2O3v2',molAl2O3v2(1))


% Calculate mole fractions
XSiO2      = molSiO2    ./ summoles;
XTiO2      = molTiO2    ./ summoles;
XAl2O3     = molAl2O3v2 ./ summoles;
XFe2O3     = molFe2O3   ./ summoles;
XFeO       = molFeO     ./ summoles + ( 2 * XFe2O3 );
XMgO       = molMgO     ./ summoles;
XCaO       = molCaO     ./ summoles;
XNa2O      = molNa2O    ./ summoles;
XK2Al2O4   = molK2O     ./ summoles;

%Check Values
%fprintf('%s %12.6f\n','summoles',summoles(3))

% Calculate activity of Anorthite in Plagioclase
%
AnTerm1    =  ( -0.000000000010696          * (TempC.^4) ) + ( 0.00000004945054       * ( TempC.^3) ) - ( 0.00008207491      * ( TempC.^2) ) + ( 0.05275448      * TempC )      - 7.914622 ;
%
AnTerm2    =  ( 0.000000000000219441        * (TempC.^5) ) - ( 0.00000000106298       * ( TempC.^4) ) + ( 0.00000200317      * ( TempC.^3) ) - ( 0.00182716      * (TempC.^2) ) + ( 0.811309      * TempC ) - 146.951 ;
% 
AnTerm3    =  ( -0.0000000000000002443527   * (TempC.^6) ) + ( 0.0000000000013185282  * ( TempC.^5) ) - ( 0.0000000028995878 * ( TempC.^4) ) + ( 0.0000033377181 * (TempC.^3) ) - ( 0.0021408319  * (TempC.^2) ) + ( 0.73454132 * TempC ) - 104.25573 ;
%
AnTerm4    =  (  0.00000000000000009642821  * (TempC.^6) ) - ( 0.0000000000005482202  * ( TempC.^5) ) + ( 0.000000001280331  * ( TempC.^4) ) - ( 0.00000157613   * (TempC.^3) ) + ( 0.001084872   * (TempC.^2) ) - ( 0.3997902  * TempC ) + 63.39784;
%
AnConstant =  ( -0.000000000000000006798238 * (TempC.^6) ) + ( 0.00000000000004008407 * ( TempC.^5) ) - ( 0.0000000000972437 * ( TempC.^4) ) + ( 0.0000001242399 * (TempC.^3) ) - ( 0.00008824071 * (TempC.^2) ) + ( 0.03310392 * TempC ) - 5.136283;
%
Aan        =  ( AnTerm1 .* (XAn.^4) )   +  ( AnTerm2 .* (XAn.^3) )  +  ( AnTerm3 .* (XAn.^2) )  + (AnTerm4 .* XAn) + AnConstant ;

%Check Values
%fprintf('%s %12.6f\n','AnTerm1',AnTerm1(1))
%fprintf('%s %12.6f\n','AnTerm2',AnTerm2(1))
%fprintf('%s %12.6f\n','AnTerm3',AnTerm3(1))
%fprintf('%s %12.6f\n','AnTerm4',AnTerm4(1))
%fprintf('%s %12.6f\n','AnConstant',AnConstant(1))
%fprintf('%s %12.6f\n','Aan',Aan(1))

% Calculate activity of Albite in Plagioclase
%
AbTerm1    =  ( 0.0000000000311852        * (TempC.^4) ) - ( 0.000000124156        * (TempC.^3) ) + ( 0.000181664        * (TempC.^2) ) - ( 0.110128        * TempC ) +   19.9603;
%
AbTerm2    =  ( -0.00000000006683324      * (TempC.^4) ) + ( 0.000000264682        * (TempC.^3) ) - ( 0.0003813719       * (TempC.^2) ) + ( 0.2230349       * TempC ) -   35.75847;
%
AbTerm3    =  ( 0.000000000045560722      * (TempC.^4) ) - (0.00000017904633       * (TempC.^3) ) + ( 0.00025203154      * (TempC.^2) ) - ( 0.13902056      * TempC ) +   17.479424;
%
AbTerm4    =  ( -0.00000000001038026      * (TempC.^4) ) + (0.00000004046137       * (TempC.^3) ) - ( 0.00005510837      * (TempC.^2) ) + ( 0.02766904      * TempC ) -    0.9465208;
%
AbConstant =  ( 0.00000000000000003664241 * (TempC.^6) ) - ( 0.0000000000002106807 * (TempC.^5) ) + ( 0.0000000004975989 * (TempC.^4) ) - ( 0.0000006178773 * (TempC.^3) ) + ( 0.0004252932 * (TempC.^2) ) - ( 0.1537074 * TempC ) + 22.74359;
%
Aab        =  ( AbTerm1 .* (XAb.^4) )  +  ( AbTerm2 .* (XAb.^3) )  +  ( AbTerm3 .* (XAb.^2) )  +  (AbTerm4 .* XAb)  +  AbConstant;

%Check Values
%fprintf('%s %12.6f\n','AbTerm1',AbTerm1(3))
%fprintf('%s %12.6f\n','AbTerm2',AbTerm2(3))
%fprintf('%s %12.6f\n','AbTerm3',AbTerm3(3))
%fprintf('%s %12.6f\n','AbTerm4',AbTerm4(3))
%fprintf('%s %12.6f\n','AbConstant',AbConstant(3))
%fprintf('%s %12.6f\n','Aab',Aab(3))
%fprintf('%s %12.6f\n','XAb',XAb(3))

% Thermodynamic data- Volumetric Data
TempK = TempC + 273;
%
VolAbxtal     = 100.570 * exp ( 0.0000268 * ( TempK - 298 ) ); %volume of albite crystal
%
VolAnxtal     = 100.610 * exp ( 0.0000141 * (TempK - 298) );   %volume of anorthite crystal
%
VliqAb    = 112.715 + ( 0.00382 * (TempK - 1373 ) ) ;      %volume of liquid at 1 bar
%
VliqAn    = 106.300 + ( 0.003708 * (TempK - 1673) ) ;      %volume of liquid at 1 bar
%
dVdPliqAb = (0.75 * (-1.843) + 0.125 * (0.685 + 0.0024 * ( TempK - 1673 ) ) + 0.125 *( -2.384 - 0.0035 * (TempK - 1673) ) ) * 4/10000; %the change in volume of albite liquid with temperature
%
dVdPliqAn = (0.50 * (-1.906) + 0.250 * (-1.665) + 0.25 * ( 0.295 - ( 0.00101 * (TempK - 1673 ) ) ) ) * 4/10000; %the change in volume of anorthite liquid with temperature
%
deltaValbite_J      = 0.1 * ( VliqAb - VolAbxtal ) .* (Pbars-1); %first part of the solution to the pressure integral for albite; AKA deltaV(P-1)albite(J) 
%
deltaValbite2_J     = 0.1 * ( ( dVdPliqAb - 0.0000167 ) / 2) .* (Pbars.^2); %second part of the solution to the pressure integral for albite; AKA dV/dP^2albite(J)
%
deltaVanorthite_J   = 0.1 * (VliqAn - VolAnxtal ).*(Pbars-1); %first part of the solution to the pressure integral for anorthite; AKA deltaV(P-1)anorth(J)
%
deltaVanorthite2_J  = 0.1 * ( ( dVdPliqAn - 0.0000116 ) / 2) .* (Pbars.^2); %second part of the solution to the pressure integral for anorthite;AKA dV/dP^2anorthtie(J)
%
intdeltaVAn         = deltaVanorthite_J  +  deltaVanorthite2_J; % the volume change between anorthite liquid and cyrstal at temperature and pressure
%
intdeltaVAb       = deltaValbite_J     + deltaValbite2_J; % the volume change between albite liquid and cyrstal at temperature and pressure
%
intdeltaV         = - deltaValbite_J   -  deltaValbite2_J + deltaVanorthite_J + deltaVanorthite2_J;% the total change in volume between anorthite and albite crystal and liquid

%Check Values
%fprintf('%s %12.6f\n','intdeltaVAn ',intdeltaVAn (3))


% Thermodynamic data- Calorimetric Data
%
XAbliq  = ( (XNa2O.^0.5) .* (XAl2O3.^0.5) .* (XSiO2.^3) ) * 18.963     ;% mol fraction of albite in the liquid
%
XAnliq  = ( XCaO .* XAl2O3 .* (XSiO2.^2) ) * 64                        ;% mol fraction of anorthite in the liquid
%
XAn_AnAb_liquid = XAnliq ./ (XAbliq + XAnliq) * 100                    ;%liquid An #
%
ln      =  @(x)(log(x))                                                ;%utilizing natural log function
%
lnK     = ( ln(Aab) - ln(Aan) ) + ( ln(XAnliq) - ln(XAbliq) )          ;%equilibrium constant
%
volterm = intdeltaV ./ (8.3144 * TempK);                                %overall change in volume considering the effects of temperature and pressure divided by the constant R and T in kelvin
%
intdeltaCp_An       = ( ( -5.77 ) * (TempK - 1830 ) ) - ( (-3734 / 0.5) * (TempK.^0.5 - 1830^0.5) ) - ( (317020000 / -2) * (TempK.^-2 - 1830^-2) ); %AKA int ?Cp (An)the integral of the heat capacity of anorthite crystal
%
deltaHfus_An        = 142406  +  intdeltaCp_An; %AKA ?H fus (An)
%
intdeltaCp_Ab       = ( ( -35.64 ) * (TempK - 1373) ) - ( ( -2415.5 / 0.5) * ( TempK.^0.5 - 1373^0.5 ) ) - ( (789280) * ( (TempK.^-1) - (1373^-1) ) ) - ( ( 1070640000 / -2 ) * ( (TempK.^-2) -(1373^-2) ) );%AKA int ?Cp (Ab); the integral of the heat capacity of albite crystal
%
deltaHfus_Ab        = 64500 + intdeltaCp_Ab; %AKA ?H fus (Ab)
%
Cap_deltaH_AnAb     = deltaHfus_An - deltaHfus_Ab; %the change in enthalpy for the exchange reaction between anorthite and albite
%
Cap_deltaH_RT       = Cap_deltaH_AnAb ./ ( 8.3144 * TempK);
%
intCp_divdedTfus_An = -5.77 * ( log( TempK ) - log( 1830 )) + ( 3734 / -0.5 ) * ( (TempK.^-0.5) - (1830^-0.5) ) + (0/-2) * ((TempK.^-2) - (1830^-2)) - ( 317020000 / -3 ) * ( (TempK.^-3) - (1830^-3) );%AKA int ?Cp/T fus (An)
%
deltaSfus_An        = 77.82 + intCp_divdedTfus_An; %?S fus (An)
%
intCp_divdedTfus_Ab = -35.64 * ( log(TempK) - log(1373) ) + ( 2415.5 / -0.5) * ( (TempK.^-0.5) - (1373^-0.5) ) + ( 7892800 / -2 ) * (TempK.^-2 - 1373^-2) - (1070640000/-3) * ( (TempK.^-3) - (1373^-3) ); % AKA int ?Cp/T fus (Ab)
%
deltaSfus_Ab        = 47 + intCp_divdedTfus_Ab; %?S fus (Ab)
%
CapDeltaS_R         = (deltaSfus_An - deltaSfus_Ab) / 8.3144; %?S/R
%
H_TdeltaS_An        = (deltaHfus_An - (TempK .* deltaSfus_An))/1000; %AKA ?H-T?S(An) in Kilojoules
%
H_TdeltaS_Ab        = (deltaHfus_Ab - (TempK .* deltaSfus_Ab))/1000; % ?H-T?S(Ab)
%
Gibbs_exchange      = H_TdeltaS_An - H_TdeltaS_Ab; %Gibbs free energy of exchange
%
Gibbs_exchange_RT   = 1000 * Gibbs_exchange ./ (8.31441 * TempK); %?G/RT
%
P_T                 = Pbars ./ TempK; %Pressure divided by temperature
%
H_RTminusS_R        = Cap_deltaH_RT - CapDeltaS_R;

%Check Values
%fprintf('%s %12.6f\n','Gibbs_exchange_RT',Gibbs_exchange_RT(3))

% Variables for H2O Calculation 
%
neglnK_V_G  = -1 * ( lnK + volterm + Gibbs_exchange_RT);
%
T_term      = 10000./TempK;
%
negXSiO2    = -1 * XSiO2;
%
negXTiO2    = -1 * XTiO2;
%
negXAl2O3   = -1 * XAl2O3;
%
negXFeO     = -1 * XFeO;
%
negXMgO     = -1 * XMgO;
%
negXCaO     = -1 * XCaO;
%
negXNa2O    = -1 * XNa2O;
%
negXK2Al2O4     = -1 * XK2Al2O4;
%
wtH2Ocalc   = - 17.3204203938587 + (0.389218669342048 * neglnK_V_G) + (2.98588693659695 * T_term) + (7.8289140199477 * negXSiO2) - (50.1063951084878 * negXAl2O3) + (14.114740308799 * negXFeO) + (23.9996276026497  * negXMgO) - (15.9003801663855 * negXCaO) + (18.6326909831708 * negXNa2O) + (24.0180473651546 * negXK2Al2O4) ;
fprintf('%s %12.6f\n','neglnK_V_G',neglnK_V_G(3));

% Calculate residual wt% H2O
Residual    = wtH2O - wtH2Ocalc;


% OUTPUT VARIABLES
% You may change the file name, which is currently 'results.csv' to
% whatever name you choose. The program automatically writes the csv file. 

fid = fopen('results.csv','w');
for i=1:size(input_data,1)
    fprintf(fid,'%6.4f, %6.4f, %6.4f, %6.4f, %6.2f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, \n', ...
                 XAn(i),XAb(i),Aan(i),Aab(i),TempC(i),Pbars(i),wtSiO2(i),wtTiO2(i),wtAl2O3(i),wtFe2O3T(i),wtFeOT(i),wtMgO(i),wtCaO(i),wtNa2O(i),wtK2O(i),wtH2Ocalc(i));
end
fclose(fid);

% Output all inputs, Aan, Aab, wtH2Ocalc values in csv format (no header)e
%output_data = [XAn,XAb,Aan,Aab,TempC,Pbars,wtSiO2,wtTiO2,wtAl2O3,wtFe2O3T,wtFeOT,wtMgO,wtCaO,wtNa2O,wtK2O,wtH2O,wtH2Ocalc,Residual];
%csvwrite('results.csv',output_data)


